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Abstract. We provide a brief overview of recent advances and outstanding issues in simulations 
of interstellar turbulence, including isothermal models for interior structure of molecular clouds 
and larger-scale multiphase models designed to simulate the formation of molecular clouds. 
We show how self-organization in highly compressible magnetized turbulence in the multiphase 
ISM can be exploited in simple numerical models to generate realistic initial conditions for star 
formation. 
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1. Introduction 

Since most of the ISM is characterized by very large Reynolds numbers, turbulent 
motions control the structure of nearly all temperature and density regimes in the inter- 
stellar gas (Elmegreen & Scalo 2004). Because of that, turbulence is often viewed as an 
organizing agent forming and shaping hierarchical cloudy structures in the diffuse ISM 
and ultimately in star-forming molecular clouds (e.g., Vazquez-Semadeni & Passot 1999). 
Nonlinear advection dominating the dynamics of such highly compressible magnetized 
multicsale self-gravitating flows makes computer simulations practically the only tool 
to study fundamental aspects of interstellar turbulence, even though effective Reynolds 
numbers in numerical models are always limited by the available computational resource 
(e.g., Kritsuk et al. 2006). 

Over the last five years three-dimensional numerical simulations fostered the devel- 
opment of theoretical concepts concerning the interstellar medium undergoing nonlin- 
ear self-interaction and self-organization in galactic disks. One can conventionally divide 
these models designed to tackle various aspects of interstellar turbulence into three differ- 
ent classes depending on the range of resolved scales and physics included: (i) mesoscale 
models that cover evolution of multiphase ISM in volumes with linear size of a few-to-ten 
kpc and resolve the flow structure down to a fraction of 1 pc (e.g., galactic fountain mod- 
els developed by de Avilles & Breitschwerdt 2005-07); (ii) sub-mesoscale models resolving 
the scale-height of the diffuse Hi (~ 100 pc) and usually limited to only warm-to-cold 
neutral phases (WNM and CNM) (e.g., Kissmann et al. 2008; Gazol et al. 2009; Gazol 
& Kim 2010; Seifried et al. 2010); and (hi) microscale models for molecular cloud (MC) 
turbulence that assume an isothermal equation of state and deal with < 10 pc-sized sub- 
volumes within MCs. Global galactic disk models (e.g., Tasker & Bryan 2006-08, Wada 
2008 and references therein) which represent the future of direct ISM turbulence mod- 
eling, are currently resolving scales down to ~ 10 pc, i.e. insufficient to properly follow 
the thermal structure of self-gravitating multiphase ISM. 

Mesoscale models of supernova-powered (SNe) galactic fountain have demonstrated the 
important role of dynamic pressure in the ISM that keeps large fractions of the gas mass 
out of thermal equilibrium and elevates gas pressures of GMCs to the observed levels 
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even without direct action of self-gravity (Korpi et al. 1999; Mac Low et al. 2005; de 
Avilles & Breitschwerdt 2005-07; Joung et al. 2006-09). They also show that the effective 
integral scale of the SN-driven turbulence (~ 75 pc) is about half the scale height of the 
Hi gas in the inner Galaxy [100 — 150 pc (Malhotra 1995)] and outline a general picture 
of probability distributions for the mass density, magnetic field strength, and thermal 
pressure in the turbulent ISM in disk-like galaxies. 

Supersonic isothermal turbulence simulations in periodic boxes representative of the 
microscale models provided many important insights into the physics of interstellar turbu- 
lence and helped to guide the interpretation of observations. These numerical experiments 
highlighted the importance of nonlinear advection as a major feature of compressible tur- 
bulence (Pouquet et al. 1991). At high Mach numbers, turbulent flows are dominated by 
shocks; therefore the velocity spectra are steeper than the Kolmogorov slope of —5/3 
and closely resemble the Burgers —2 scaling (Kritsuk et al. 2006b). The physics of three- 
dimensional supersonic turbulence is, however, quite different from burgulence (Frisch 
& Bee 2001) as the solenoidal velocity component always remains dominant (Pouquet 
et al. 1991; Pavlovski et al. 2006; Kritsuk et al. 2007; Pan et al. 2009; Schmidt et al. 
2009; Kritsuk et al. 2010). At sonic Mach numbers M s > 3, strong shock interactions 
and associated nonlinear instabilities create sophisticated multiscale pattern of nested 
U-shaped structures in dynamically active regions morphologically similar to what is 
observed in molecular clouds (Kritsuk et al. 2006a). Scaling of the first-order velocity 
structure functions Si(5u) ~ £°- 5i (where Su(£) = u(x) - u(x + U), S p (6u) = ([<5u(£)F) 
and (. . .) indicates averaging over an ensemble of random point pairs separated by the lag 
£) obtained in simulations (Kritsuk et al. 2007) is similar to the velocity scaling observed 
in molecular clouds Si(Su) ~ £ ' 56 (Heyer & Brunt 2004). Simulations also support the 
concept of (lossy) energy cascade in compressible turbulence (e.g., Vazqucz-Semadeni et 
al. 2003), suggesting that the kinetic energy directly lost in shocks constitutes a small 
fraction of the total energy dissipation. The fact that the Richardson-Kolmogorov cascade 
picture does approximately hold for supersonic turbulence follows from the linear scaling 
of the third-order structure function of the mass- weighted velocity, S^Sf/pu) ~ £, indi- 
cating constant turbulent energy transfer rate across the hierarchy of scales (Kritsuk et 
al. 2007; Kowal & Lazarian 2007; Schwarz et al. 2010). The power spectra of tfpu, accord- 
ingly, demonstrate the Kolmogorov scaling independent of the Mach number (Kritsuk et 
al. 2007; Schmidt et al. 2008; Kritsuk et al. 2009; Federrath et al. 2010; Price & Feder- 
rath 2010). It seems that this result can be also extended to supersonic MHD turbulence, 
where the incompressible 4/3-law of Politano & Pouquet (1998) also approximately holds 

in its scaling part, Sf 3 = (SZJ (l)[SZf (£)] 2 ^ ~ £ (here SZ {1 (£) = [Z(x + U) - Z(x)] • e, 

and e£ is the displacement vector) , if reformulated in terms of the mass- weighted Elsasser 
fields Z^ = p 1 / 3 (u± TS/ ^JAirp) (Kritsuk et al. 2009). The presence of magnetic field effec- 
tively reduces compressibility of the gas making the velocity spectra more shallow with 
slopes approaching the Iroshnikov-Kraichnan index of —1.5 in trans- Alfvcnic flows. The 
observed spectral slope of about —1.8 for the Perseus molecular cloud is thus consistent 
with the super- Alfvenic turbulence regime dominant in that cloud (Padoan et al. 2006). 

Recent results from numerical experiments on highly compressible turbulence stim- 
ulated theorists to reconsider the steady-state statistics of turbulence in the inertial 
interval. Falkovich et al. (2010) have shown that the Kolmogorov 4/5-law is a particular 
case of the general relation on the current-density correlation function. They derived an 
analog of the flux relation for compressible turbulence that can be used as a test for 
direct numerical simulations and as a guide for the development of subgrid scale models 



for astrophysical turbulence (Schmidt & Federrath 2010) 
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Most of the recent sub-mesoscale models belong to a class of so-called converging (or 
colliding) flows of diffuse Hi originally developed to study thermal, dynamic, and gravita- 
tional instabilities in shock-bounded slabs (Hunter et al. 1986; Vishniac 1994; Walder & 
Folini 1998; Folini et al. 2010). These models remain popular as a framework to directly 
simulate star formation in molecular clouds (Vazquez-Semadeni et al. 2006; Hennebelle 
& Inutsuka 2006; Vazquez-Semadeni et al. 2007; Hennebelle et al. 2008; Inoue & Inutsuka 
2008-09; Heitsch et al. 2008-09; Banerjee et al. 2009; Niklaus et al. 2009; Audit & Hen- 
nebelle 2010; Rosas-Guevara et al. 2010). Recent numerical experiments with converging 
flows have demonstrated strong sensitivity of results to adopted initial and boundary con- 
ditions as well as to model parameters that control the density of colliding gas streams, 
mean thermal pressure, orientation and strength of the mean magnetic field, levels and 
character of "turbulence" at infinity, etc. All these parameters live their unique imprints 
in the statistics of derived stellar populations and any comprehensive parameter study 
based on computational modeling in this framework would be prohibitively expensive. 

One way to circumvent these difficulties is to exploit Prigogine's concept of self- 
organization in nonequilibrium nonlinear dissipative systems (Nicolis & Prigogine 1977) 
in application to the ISM (e.g., Biglari & Diamond 1989). With this approach, one can 
use interstellar turbulence as an agent that imposes "order" in the form of coherent 
structures and correlations between various flow fields emerging in a simple periodic box 
simulation when a statistical steady state develops. In this case, the initial conditions 
are no longer important, instead the steady state would provide the "correct" turbulent 
initial conditions for star formation when self-gravity is turned on. While this idea is not 
new0 it remained largely undeveloped so far. In the following sections we will discuss 
this concept in more detail and report first results from a series of MHD simulations of 
turbulent multiphase ISM with the piecewise parabolic method on a local stencil (PPML; 
Ustyugov et al. 2009). 



2. Self-organization in the magnetized multiphase ISM 

In out numerical experiments, we treat the ISM as a turbulent, driven system, with 
kinetic energy being injected at the largest scales by supernova explosions, shear associ- 
ated with differential rotation of the galactic disk, gas accretion onto the disk, etc. (Mac 
Low & Klessen 2004; Klessen & Hennebelle 2010). This kinetic energy is then being 
transferred from large to small scales in a cascade-like fashion. As our models include a 
mean magnetic field, Bq, some part of this kinetic energy gets stored in the turbulent 
magnetic field component, b, generated by stretching, twisting, and folding of magnetic 
field lines. The ISM is also exposed to the far-ultraviolet (FUV) background radiation 
due to OB associations of quickly evolving massive stars that form in molecular clouds. 
This FUV radiation is the main source of energy input for the neutral gas phases and this 
volumetric thermal energy source is in turn balanced by radiative cooling (Wolfire et al. 
2003). The ISM is thus exposed to various energy fluxes, and self-organization arises as a 
result of the relaxation through nonlinear interactions of different physical constituents 
of the system subject to usual MHD constraints in the form of conservation laws. In this 
picture, molecular clouds with their hierarchical internal structure form as dissipative 
structures that represent active regions of highly intermittent turbulent cascade that 
drain the kinetic energy supplied by the driving forces. 



f See, for instance, summary of the panel discuss ion on Phases of the ISM during the 1986 
Grand Teton Symposium in Wyoming (Shull 1987) 
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Table 1. Model parameters. 
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Figure 1. Three snapshots of projected gas density in model A taken at t = 2, 3, and 4 Myr. 
The white-blue-yellow colors correspond to low-intermediate-high projected density values. 



3. Modeling the formation of molecular clouds 

To illustrate these ideas, we consider a set of simple periodic box models, which ignore 
gas stratification and differential rotation in the disk and employ an artificial large-scale 
solenoidal force to mimic the supply of kinetic energy from various galactic sources. This 
naturally leads to an upper bound on the box size, L, which determines our choice of 
L = 200 pc. Our models are, thus, fully defined by the following three parameters: the 
mean gas density in the box, no; the rms velocity, u rms .o; an d the mean magnetic field 
strength, Bq. All three would ultimately depend on L; Table 1 provides the summary of 
parameters for models A, B, C, D, and E assuming L — 200 pc. The table also gives the 
grid resolution, N, the initial values for plasma beta, /3th,o = 8ttpq/Bq, turbulent beta, 
/?turb,o = 87r Po u ?ms .o / ; an d Alfvenic Mach number, Ma.o = (47rpo) 1 ^ 2u rms,o/-So! where 
Po is the initial thermal pressure of the gas (see the phase diagram in Fig. [2] for more 
detail). 

We initiate our numerical experiments with a uniform gas distribution in the compu- 
tational domain. An addition of small random isobaric density perturbations at t = 
triggers a phase transition in the thermally bi-stable gas that quickly turns ~ 25 — 65% 
of the gas mass into the stable cold phase (CNM with temperature below T = 184 K), 
while the rest of the mass is shared between the unstable and stable warm gas (WNM). 
In models A, B, and C, CNM and WNM each contain roughly ~ 50% of the total Hi 
mass in agreement with observations (Heiles & Crutcher 2005). We then turn on the 
forcing and after a few large-eddy turn-over times the simulation approaches a statistical 
steady state. If we replace this two-stage initiation process with a one-stage procedure 
by turning the driving on at t — 0, the properties of the steady state remain unchanged. 
Figure Q] illustrates this evolutionary sequence for the two-stage case with three snap- 
shots of projected gas density for model A. The left panel shows two-phase medium at 
t = 2 Myr right before we turn on the forcing; the panel in the middle illustrates an early 
stage of turbulization with transient "colliding flows" at t = 3 Myr. The right panel 
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Figure 2. Time evolution of the rms magnetic field strength, _B rms , and its turbulent component, 
ferms, for models A, B, C, and D (left panel). Phase diagram (thermal pressure versus gas density) 
for model B at t = 5 Myr (right panel). 




Figure 3. Time-average density distributions for fully developed turbulence in models A, B, 
and C (left panel) and time-average mass- weighted thermal pressure distributions for models A, 
B, C, and D (right panel); see text for more detail. 



shows the projected density at t = 4 Myr for a statistically developed turbulent state. 
Molecular clouds can be seen in the right panel as filamentary brown-to-yellow structures 
(note that these are morphologically quite different from the transient dense structures 
in the middle panel). The rms magnetic field is amplified by the forcing and saturates 
when the relaxation in the system results in a steady state, see Fig. [2] The level of satu- 
ration depends on Bq and on the rate of kinetic energy injection by the large-scale force, 
which is in turn determined by u rms and uq. This level can be easily controlled with the 
model parameters. In the saturated regime, models A and B tend to establish energy 
equipartition (i?K ~ Em), while the saturation level of magnetic energy in model C is 
a factor of ~ 3 lower than the equipartition level. The mean thermal energy also gets a 
slight boost due to forcing, but remains subdominant in all the models. A typical phase 
diagram is shown in Fig. [2l The contours indicate constant levels of volume fraction for 
different regimes of the thermal pressure, pth, and density, n, separated by factors of 2. 
About 23% of the domain volume is filled with the stable warm phase at T > 5250 K, 
the stable cold phase (T < 184 K) occupies ~ 7%, and ~ 70% of the volume resides 
in the thermally unstable regime at intermediate temperatures. The big orange dot at 
the center indicates the (forgotten) initial conditions for models A, B, C, and E. The 
phase diagram indicates that turbulence supports an enormously wide range of thermal 
pressures and also that pth in the molecular gas (n > 100 cm -3 ) is higher than that in 
the diffuse ISM, even though self-gravity is ignored in the model. 
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Figure 5. Distributions of the AlfVenic Mach number and the cosine of the alignment angle 
vs. density for a snapshot from model A at t = 5 Myr. 



Figure [3] (left) shows the time-average density PDFs for models A, B, and C in the 
steady state. The effect of magnetic field on the density PDF is apparently very weak on 
average, and the high-density part of the PDF can be well approximated by a lognormal 
as in the isothermal case. The signature of the two stable thermal phases in the PDF is 
smeared by the (relatively) high level of turbulence, u rm s(100 pc) = 16 km/s (Brunt & 
Heyer 2004) , but the overall shape of the distribution is not lognormal. The distribution 
of thermal pressure for models A, B, and C spans about 6 dex leaving no room for 
the old pressure-supported cloud picture in the violent ISM. All distributions match the 
characteristic pressure typical for the Milky Way disk at the solar radius and show only 
weak dependence on Bq , while the width of the distribution remains sensitive to u rms and 
uq. Figure [3] ( right) shows how the mass- weighted pressure distributions obtained in our 
models compare with the distribution reconstructed from high-resolution UV spectra of 



hot stars in the HST archive (Jenkins & Tripp 2010) It seems that models with zt rms = 
16 km/s reproduce both the shape and the width of the observed distribution quite nicely, 
while a lower turbulence level in model E makes that model distribution too narrow. 

These numerical experiments allow us to probe the levels of magnetic field strength in 
molecular clouds that form self-consistently in the magnetized turbulent diffuse ISM. In 
the left panel of Fig.|4j we show a scatter plot of magnetic vs. thermal pressure for model B 
at t = 5 Myr. The black contour lines show the distribution for the whole domain, which is 
centered at (3 t h ~ 0.1. The subset of cells representing the molecular gas (T < 100 K and 
n > 100 cm" 3 , color contour plot) shows very similar mean values of /3th- This indicates a 
target plasma beta for realistic isothermal molecular cloud turbulence simulations. The 
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right panel of the same Figure shows a scatter plot of magnetic vs. dynamic pressure 
for the same snapshot from the same model. The distribution for the whole domain is 
centered at /3turb ~ 1 because kinetic and magnetic energy levels are close to equipartition 
on average. At the same time, the subset of cells representing the molecular gas (color 
contours) shows a distribution centered at /3turb ~ 30 meaning that turbulence in the 
molecular gas is super- Alfvenic. Figure [5j left panel, shows the distribution of Alfvenic 
Mach number, Ma, as a function of density for the strongly magnetized model A that 
further supports this result. There is a clear positive correlation, Ma ~ n 0A , indicated 
by the dashed line and most of the dense material (n > 100 cm -3 ) clearly falls into the 
super- Alfvenic part of the distribution. 

A key to understanding the origin of this super- Alfvenic regime in the cold and dense 
molecular gas lies in the process of self-organization in magnetized ISM turbulence that 
we briefly introduced in Section [2] The statistical steady state that our models attain 
on a time-scale of a few million years is characterized by a certain degree of alignment 
between the velocity and magnetic field lines. Figure[5j right panel, shows the distribution 
of the cosine of the alignment angle, cos 9 = B ■ u/(Bu), for model A at t = 5 Myr. The 
contours indicate a saddle-like structure of this probability distribution with a strong 
alignment regime (cos 9 = ±1) in the WNM at and around n ~ 1 cm -3 . This means that 
compressions in the WNM gas, which is on average trans-Alfvenic (e.g., Ma £ [0.6,0.9] 
in model A), occur preferentially along the field lines. If molecular clouds form in the 
turbulent ISM via large-scale compression of the diffuse Hi, then turbulence in such 



molecular clouds can only be super- Alfvenic, see also Padoan et al. (2010) 



4. Conclusion 

Rapid development of computational astrophysics in the recent years enabled progress 
in understanding the basics of interstellar turbulence. These new advances will help us 
to move forward with direct star formation simulations from turbulent initial conditions. 
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